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Introduction 

The purpose of this research is to develop a rigorous theory and corresponding com- 
putational algorithms for a variety of problems regarding the analysis of composite beams 
and plates. The modeling approach is intended to be applicable to both static and dy- 
namic analysis of generally anisotropic, nonhomogeneous beams and plates. The major 
part of the effort during this first reporting period has been devoted to development of a 
theory for analysis of the local deformation of plates. In addition, some work has been 
performed on global deformation of beams. Because of the strong parallel between beams 
and plates, we will treat the two together as thin bodies, especially where we believe it 
will aid in clarifying to the reader the meaning of certain terminology and the motivation 
behind certain mathematical operations. 


Background 

The static and dynamic analysis of beams and plates is of fundamental importance in 
many engineering problems. In design and analysis of modem aerospace systems, subsys- 
tems which consist of laminated plates or composite beams are often encountered such as 
wing, fuselage, and rotor blade structures. Classical beam and plate theories are known 
to be adequate for many applications. For beams or plates with anisotropic and nonho- 
mogeneous construction, however, these theories suffer from several sources of inaccuracy. 
For example, composite beams — even those which are slender — can be very flexible in 
shear. When such beams are analyzed by classical theories, which generally ignore shear 
deformation, certain kinematical quantities and associated constitutive couplings are ab- 
sent which are known to be important. Rehfield, Atilgan, and Hodges (1990) concluded 
that such phenomena can influence the global deformation in thin-walled beams designed 
for extension-twist coupling by as much as a factor of two! One of the chief mechanisms for 
this large effect is a coupling between bending and shear deformation due to anisotropic 
materials (see, Rehfield and Atilgan, 1989). 

This conclusion also holds for beams which are not necessarily thin-walled; how- 
ever, the cross sectional analysis of general nonhomogeneous, anisotropic beams cannot 
be treated without finite elements. Unlike classical beam or plate theories, theories for 
composite beams or plates must contain some means of calculating the elastic properties. 
These properties are not merely material moduli multiplied by certain integrals over the 
section of the beam. A review by Hodges (1990a) covers most of the literature dealing 
with modeling of beams prior to 1988. A similar treatment of composite plates, in which 
the determination of properties is carried out as a separate analysis, appears to be missing 
from the literature. 

The analysis of deformation for composite thin bodies normally must be done as 
an iterative process. There are local deformation analyses, which determine the elastic 
constants and details of local deformation in terms of global deformation parameters. Also 
there are global analyses, which determine the global (beam- or plate-like deformation such 
as bending, extension, torsion, and shear) deformation. One could proceed as follows. (1) 
calculate the local deformation and elastic constants for the undeformed structure; (2) use 
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the resulting elastic constants in a global deformation analysis; (3) recalculate the local 
deformation for the deformed structure; (4) stop the process if the elastic constants do 
not change more than some tolerance or else go back to step 2. Most plates and many 
beams do not exhibit a sufficiently large local deformation (which we will call warping) to 
require more than one calculation of properties. In other words, elastic constants which 
are determined based on the undeformed state may be sufficient even for geometrically 
nonlinear analyses. 

There is another aspect of research which we originally proposed, that being the use of 
space-time finite elements to treat the global dynamics of beams. In the time between the 
submittal of our proposal to the Aerostructures Directorate and its being funded, another 
of our proposals was funded. It had been under review at the National Science Foundation 
for a long time; among other things, it concerns the dynamics of beams by space-time finite 
elements. Since both of these grants have this common subject matter as subsets of their 
programs, we intend to treat this portion of the research as jointly funded. Prof. David 
A. Peters and Dr. Weiyu Zhou have contributed to the NSF work, and thus, indirectly to 
this project. 

Most finite element procedures for time-dependent phenomena are based on semi- 
discretizations; finite elements are used in space to reduce to a system of ordinary differen- 
tial equations. This kind of procedure is widely used in practice and fairly well understood. 
Space-time finite elements have been rarely used in solutions of engineering dynamic prob- 
lems. Classic time integration methods are usually included in computational procedures. 
Recent developments of the space-time finite element method allow application of approx- 
imation techniques to the spatial and temporal domains. Special schemes lead to highly 
efficient algorithms that reduce both memory requirements and number of arithmetical 
operations. 

The use of space-time finite elements presents yet another duality. The static and 
frequency-domain global deformation analyses for plates are solved on the 2-D domain of 
the plate; the space-time dynamics of beams can be cast as a 2-D problem with simulta- 
neous spatial and temporal discretizations. 

Previous Work Related to Beams 

Although there has been much work published in the literature on beams and plates, 
a unified approach such as we are undertaking has not received very much attention. 
Furthermore, most of what does exist in the literature applies to beams only. Thus, before 
getting into plate references, we will describe the analogous type of analysis for beams with 
the hope of clarifying what we have begun to do for plates. 

Berdichevsky (1981) was the first to prove, from variational-asymptotic considera- 
tions, that nonlinear analysis of beams can be split into two separate problems: the local 
deformation is a linear 2-D problem, and the global deformation is a nonlinear 1-D prob- 
lem. The global deformation analysis of beams can be undertaken by use of 6 displacement 
variables associated with a cross section of a beam. The variables describe displacement 
of a point in the cross section and rotation of the cross section as a rigid body. There are 
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also 6 generalized strains (extension and two shearing strains at the reference line, twist, 
and two bending “curvatures”) and 6 stress resultants (axial force and two components 
of shear force, twisting moment, and two components of bending moment); these must 
be related by elastic constants. This way of describing the global deformation requires at 
most 21 elastic constants (a symmetric 6x6 matrix). 


Beam Local Defojmwtion To find these constants, all possible local deformations (t.e., 
warping) of the cross section must be taken into account, although they can be assumed 
small. Here warping refers not only to out-of-plane distortion of the cross section during 
torsional deformation as in classical theories, but in-plane and out-of-plane deformation; 
these are fully coupled when the beam is nonhomogeneous and anisotropic. 

As pointed out by Hodges (1990a), the work by Giavotto et al. (1983) is the most gen- 
eral of all the published works. Giavotto et al. (1983) developed a finite element approach 
for determining the elastic constants for an arbitrarily nonhomogeneous, anisotropic beam. 
This approach makes use of a two-dimensional finite element mesh representative of the 
beam cross section geometry and material properties. The results, in addition to the ma- 
trix of elastic constants, include the distribution of warping displacements per unit values 
of each of the stress resultants (section forces and moments in the cross section basis) and 
the three-dimensional stress and strain values throughout the cross section per unit values 
of the stress resultants. 

Because of the generality of the work of Giavotto et al. (1983), we have adapted a ver- 
sion of the code developed by Giavotto, Borri, and their associates for obtaining the elastic 
constants for anisotropic beams. Its only shortcoming is its rather long computer times for 
calculating properties of realistic composite beams. For this reason, under Army Research 
Office sponsorship, we are also developing a theory based on the variational-asymptotic 
method as formulated and applied to nonlinear analysis of isotropic shells by Berdichevsky 
(1978, 1979). This method has some promise of yielding a more computationally efficient 
algorithm for extracting the properties. 


Beam Global DefwmuUwn These elastic constants can then be used to find either linear 
or nonlinear global deformation, free-vibration modes and frequencies (see, Hodges and 
Atilgan et al ., 1989, 1990), and buckling behavior (see, Rehfield and Atilgan, 1989). Vari- 
ous modal, direct numerical integration, and finite element methods exist for this purpose. 
Because of their computational efficiency and modeling flexibility, finite element methods 
are quite popular. Displacement finite element methods for geometrically nonlinear behav- 
ior of beams, however, require numerical quadrature of highly nonlinear functions of the 
beam deformation. This tends to make the numerical solution procedure quite inefficient. 
On the other hand, with a mixed method such as that of Hodges (1990b), numerical ele- 
ment quadrature can be avoided (as long as applied load terms that are explicit in the axial 
coordinate are integrable in closed form). Such a nonlinear global analysis gives the beam 
displacements, rotations, extensional strain, shear strains, twist, bending curvatures, and 
sectional forces and moments to a comparable level of accuracy. It should be noted that 


4 


one can then use these results for forces and moments to find pointwise stress or strain 
levels throughout the cross section using the cross sectional finite element mesh. 

Beam behavior in the time domain can be calculated by finite elements. For the 
dynamics of beams, recent work in the field of space-time finite elements applied to struc- 
tural dynamics are described and reviewed by Bajer and Bonthoux (1988). There have 
been several developments in this area for time domain analysis of simple linear oscillators 
and rigid bladed helicopter rotors; for example Borri et al. (1985, 1988), Borri (1986), and 
Peters and Izadpanah (1988). The results indicate that finite elements in time provide a 
way of determining the dynamic behavior of a deformable body undergoing time- dependent 
loading. We have concentrated our work on the determination of the dynamic response 
of flexible structures by simultaneous discretization of the spatial and temporal domains. 
The main purpose is to determine if this methodology can be made feasible. 


Previous ’Work Related to Plates 

Reissner (1985) departed from a three-dimensional statement of the problem recalling 
the fact that in the absence of the adjective thin , plate theory would be no more than 
a class of boundary value problems in three-dimensional elasticity. He gave an excellent 
sociological-historical survey with his interpretations of the nature of approximations and 
their consequences. He also touched upon the qualitative differences in the modeling of 
laminated plates. Since we believe that a consistent plate theory should include all possible 
deformations, our starting point is, naturally, the three-dimensional kinematics. 

The analysis of laminated plates has attracted an enormous amount of attention. A 
description of the plate problem, analogous to the beam, has not been published to the 
best of our knowledge. Although the review articles by Bert (1984) and Noor and Burton 
(1989) contain hundreds of references which treat the laminated plate problem, all of them 
use some sort of explicit, analytic, through-the-thickness assumptions for displacement. 
While this is not incorrect, simple, low-order theories of this type will not yield correct 
stresses and strains through the thickness. The reason for this is rather obvious: in a 
laminated plate where each layer has distinct material properties, the actual displacement 
is not and cannot be analytic. When the assumed displacement is differentiated to yield 
the strain, the approximate strain will also be analytic while the actual strain may be 
discontinuous. Higher-order theories will be more accurate, but at a potentially high cost. 

There axe alternatives. One is to use a family of approaches as outlined and reviewed 
by Librescu and Reddy (1989) in which breaking the plate into finite elements through the 
thickness is advocated. This will yield the correct answer, but, again, at a potentially high 
computational cost. Another alternative is to derive a local deformation theory similar to 
that for the beam described above. 

Present Approach 

Plate Modeling In this research, we are developing such a computational method for 
determining the elastic constants for laminated plates. The approach is very similar to 
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that described above for beams. Yet, to the best of our knowledge, the present approach 
has never been attempted. We originally set out to solve the interior St.-Venant solution 
for the plate by finite elements in a manner that is strictly analogous to the approach of 
Giavotto et al. (1983) for the beam problem. After applying the finite element method, 
a set of linear equations were obtained that were very similar to those of Giavotto et 
al. (1983). Solution of these equations in a manner similar to the earlier ones, however, 
turned out to be difficult. Thus, we tinned to an approximate solution based on the 
variational-asymptotic method. 

The domain of the local deformation problem for the beam is planar (2-D), just as the 
global deformation problem for the plate is. On the other hand, the global deformation 
problem for the beam is solved along a line (1-D), just as the local deformation problem for 
the plate is. For the local deformation of the plate, instead of an arbitrary interior cross 
section as with the beam, we work with an arbitrary interior normal line element of the 
plate (a line of material points normal to the reference surface of the plate). The tractions 
acting on this fine element are used to obtain a variational principle governing the local 
stress resultants and deformation of the fine element. The variational principle leads to a 
symmetric 8x8 matrix of elastic constants (for a total of 36) based on the linear relation 
between the 8 stress resultants and 8 generalized strains. 

As with the beam problem, the elastic constants will be determined from a finite 
element code that is linear. This code will enable us to calculate the elastic constants for 
an arbitrary laminated plate. In this report, after presentation of the theory, we present 
some preliminary results. These results, which essentially duplicate classical theory, were 
obtained to check out the methodology and the code. 


Dummies of Beams An important step towards obtaining a general and consistent form 
of beam elastodynamic equations was taken by Hodges (1990b). Therein, geometrically 
nonlinear beam elastodynamic equations as derived from Hamilton’s Principle; also, using 
appropriate Lagrange multipliers a mixed variational formulation suitable for space-time 
finite element discretization was developed. In order to exploit the usefulness of space-time 
finite elements we decided to start from the very basic linear equations for longitudinal 
dynamics of a beam, a special case of Hodges (1990b), and their solutions. In this report, 
after a brief treatment of the theory, there are a few preliminary results presented. 

The report closes with a description of work to be undertaken during the next reporting 
period. 

Unified Variational Formulation for Anisotropic Plates 

Our starting point is the nonlinear kinematics of deformation for plates. We will 
develop the three-dimensional Biot strain field based upon the kinematical development of 
Danielson (1989). After this, we will formulate the strain energy. Finally, approximations 
of this strain energy function will be discussed and asymptotically correct solutions for 
the through- the-thickness analysis of isotropic plates, as well as the corresponding elastic 
constants, will be given. 
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Kinematics 

A plate is a flexible body in which matter is distributed about a planar surface so 
that one dimension is significantly smaller than the other two. (Much of our analysis 
can be easily extended to treat shells, but herein we will consider only plates.) The 
reference surface is an arbitrary planar surface, not necessarily the mid-surface of the 
plate. Throughout the analysis, Greek indices assume values 1 or 2, Latin indices assume 
values 1, 2, and 3 and repeated indices are summed over their ranges. Let us establish 
a Cartesian coordinate system x, so that x Q denote lengths along orthogonal lines in the 
reference surface and X3 is the distance the normal to the reference surface. Let b, denote 
an orthogonal reference triad along the undeformed coordinate lines. The position vector 
to an arbitrary point along the normal line is 


r*(xi,x 2 ,x 3 ) = r(xi,x 2 ) + X3b 3 = xyb, (1) 

Covariant and contravariant undeformed base vectors are defined as, respectively, 


gi 


dr* 

dx{ 

1 


g =; 


--Cijk- 


dr* dr* 


( 2 ) 


2 i/g~' J '‘dxj dx k 

where g = det(g, • gy). For this analysis, both reduce to 


g, = g' = b, (3) 

In a similar manner, consider the deformed state configuration. The particle which 
had position vector r*(xi,x 2 ,x 3 ) in the undeformed plate now has position vector 
R*(xi,x 2 ,x 3 ) relative to the same point, which can be represented by 


R*(xi,x 2 ,x 3 ) = R(xi,x 2 ) + x 3 B 3 (xi,x 2 ) + ttfi(xi,X 2 ,x 3 )B;(xi,x 2 ) (4) 

where R(xi,x 2 ) = r(xi,x 2 ) + u(xj,x 2 ) and u = ti;b, is the displacement vector of the 
points on the reference surface and u>,(xi , x 2 , x 3 ) is the general local (t.e., warping) dis- 
placement of an arbitrary point on the normal line, consisting of both in- and out-of-plane 
components, so that all possible deformations are considered (Figs. 1, 2). The relationship 
between B* and bi is given by 


B,(xi , x 2 ) — Ciy(xi, x 2 )by 


( 5 ) 


where C(xi , x 2 ) is the matrix of direction cosines. Covariant deformed base vectors 


G,= 


dR* 

dxi 


( 6 ) 
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can be obtained by standard means. It should be noted, however, the measure numbers Wi 
provide redundant information since the normal line undergoes rigid-body displacement 
due to U{ and rigid-body rotation due to C. Therefore, some means of removing this 
redundancy must be introduced. For a plate, we can choose the unit vectors B, so that 
wi is small, at least in some sense. Setting an appropriate number of weighted average 
displacements to zero is one way to remove the redundancy. This can be conveniently done 
by the finite element method and will be dealt with below. 

Here we restrict ourselves to the case when strain and local rotation are small so that 
the three-dimensional Biot strain can be expressed as 


r-±±£-i 


( 7 ) 


where A,j — B, • G*g fc • by is the deformation gradient matrix and / is the 3 x 3 identity 
matrix. Here T* is a 3 x 3 symmetric matrix. Introducing the column matrix w with 
components iw,-, one can expressed the three-dimensional strain field as a 6 x 1 column 
matrix r = [Tn 2fi2 1^22 2Fi3 21^3 so that 


T — Tie + I 3 w t3 + Jitui -f I 2 w t2 


( 8 ) 


where ( ) t denotes the partial differentiation with respect to x and 
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and e = \j /cJ T , the intrinsic strain measure which is function of only x a . Also 7 = 
[711 2712 722 2713 2 7 23j T and k — [«u 2k\ 2 k 22 \ T are the so-called force and moment 
strain measures. Note that 711 and 722 are the extensional strains of the reference surface, 
712 is the shear strain in the plane of the reference surface, 7^3 are the transverse shear 
strains of the normal line element, Kn and k 22 axe the elastic components of the bending 
curvature, and k \ 2 is the elastic twist. The force and moment strain measures are so 
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designated because they are conjugate to the actual running stress and moment resultants, 
respectively. 

Since warping displacements are supposed to be quite small, the few nonlinear terms 
in the strain field, which couple w and e, have been neglected in Eq. ( 8 ). The form of the 
strain field is of great importance because it is now linear in 7 , k, and w and its derivatives. 
If the top and the bottom of the line element through the thickness of the plate are free 
of tractions, application of the principle of virtual work to an infinitesimal line element, 
would lead to a system of linear equations over the one-dimensional line element governing 
w. The warping could then be determined in terms of these intrinsic strain measures or 
stress resultants as in Giavotto et al. (1983). This would lead to a unique two-dimensional 
strain energy function t/(7)*0* The elastic law could then be put in a form 



( 10 ) 


where F and M axe column matrices F — [i'll -F12 ^22 ^13 ^23 \ T and M = 
|Mh M 12 M 22 J T Here F u and F 22 are the in-plane stretching stress resultants, F 12 
is the in-plane shear stress resultant, F a 3, are the transverse shear stress resultants, Mu 
and M 22 are the bending moment resultants and M\ 2 is the twisting moment resultant, 
all expressed in the Bj basis. The elastic stiffness matrix relating F and M to *y and k is 
8 x 8 . The matrix A is 5 x 5, D is 3 x 3, and B is 5 x 3. 

Even though this methodology is successful for beam formulations (for linear analysis 
see Giavotto et a/., 1983, and for nonlinear analysis see Atilgan and Hodges, 1990), it has 
not been completely resolved for plate analysis. (This analysis is outlined in the Appendix 
as far as we have been able to take it.) Therefore, we have changed our methodology from 
direct to asymptotical analysis. 


Strain Energy and Approximations 

One of the most consistent ways to obtain the constitutive law for thin-body (t.e. beam 
and plate /shell) analysis is the use of asymptotical analysis. Literature for the successful 
analysis of asymptotical methods for beams can be found in Hodges (198 i ) and Atilgan and 
Hodges (1990) and for plates in Noor and Burton (1989). In addition to the direct asymp- 
totical analysis, which is applied to differential equations , Berdichevsky (1978, 1979) devel- 
oped the variational- asymptotical analysis, which is applied to functionals. Berdichevsky 
and his co-workers applied this method successfully to beams and shells for static and 
dynamic analysis. An outline and a simple application of this method to beam analysis 
can be found in Berdichevsky (1980). In what follows we will apply this method to non- 
homogeneous and anisotropic plates to obtain the most consistent approximations of the 
constitutive law. 


Three-Dimensional Strain Energy for Anisotropic Plates 

The three-dimensional strain energy for an anisotropic plate can be written as 
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u = - I f r T Drdx 3 dA 

2 Ja Jh 


( 11 ) 


where D is the 6x6 symmetric material stiffness matrix which relates the three-dimensional 
Biot strain T to the three-dimensional Jaumann stress Z. The form of this matrix can be 
found in Jones (1975) for all possible type of material structures ranging from transversely 
isotropic case to most general anisotropy. The Jaumann stress is also arranged in a 6 x 1 
column matrix form Z =[ Zj i Z \ 2 Z 22 Z 23 Z 33 J so that 


Z = DT 


( 12 ) 


Since waxping is a three-dimensional function of all coordinates, for most general 
configurations, it is not possible to deal with it only through the thickness. Therefore, one 
should discretize the warping as follows 


u>(xi,x 2 ,x 3 ) = N (x 3 ) ^( 11 , 0 : 2 ) (13) 

Then, using this together with Eq. (8) in Eq. (11), one can obtain the strain energy as 
follows 


u =\LV m+ 

e T R T W + W T Re + W T EW+ 

e T LlW >a + W T a L Q e + W T a C^W + W T C Q W, a + 


W T a M aP W, p 


dA 


(14) 


where 


A= f H T DHdx 3 R= [ N ,T ljDHdx 3 E = f N' T ij DI 3 N' dx 3 

Jh Jh Jh 

L a = [ N T llDHdx 3 C a = [ N lT lJ DI Q Ndx 3 M af) = / N T I Q DI 0 Ndx 3 

Jh Jh Jh 


and where ( )* denotes differentiation with respect to x$. Because the description of the 
displacement is 5 times redundant, the rigid-body portion of the warping degrees of freedom 
must be removed in forming these equations. After this, the matrix E will be positive 
definite. The rigid-body portion of w can be removed by constraining the finite element 
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nodes. In order to do this, we set tt >3 = 0 at the reference surface and w a = 0 at the upper 
and lower surfaces. 

Development of Finite Element Approximations for Anisotropic Plates 

The zeroth approximation of the strain energy can be obtained by neglecting the 
warping completely. Then, the energy is only due to the global measures of the deformation 
and the stiffness matrix of the plate is just the matrix A. It can be shown that the matrix 
A gives an upper bound for the stiffnesses. 

A first approximation of this functional can be obtained by taking only the first 
two lines of the strain energy functional. (This approximation cannot be called the first 
approximation un til it is proven. We address this below.) These four terms are considered 
to be the most dominant terms in the strain energy functional. The reason for this is 
that differentiation with respect to the in-plane coordinates xi and xi will always result 
in smaller magnitudes than differentiation with respect to the thickness coordinate X 3 . 
Therefore, the energy obtained by the remainder terms should be smaller than the first 
four terms. The first approximation of the strain energy functional then reads 

U* =- f (e T Ae + e T R T W + W T Re + W T EW)dA (16) 

2 J a 

Since at the beginning we consider waxping to be an arbitrary quantity, independent of the 
globed strain measures, the Euler-Lagrange equation associated with W will be obtained 
from this functional by taking a straightforward variations of with respect to W yielding 

Re + EW = 0 (17) 


From this we obtain a relationship between our global strain measure e and the warping 
W to be 


W = -E~ 1 Re 


( 18 ) 


In order to prove that this is the first approximation we need to find the second approxi- 
mation. It can be shown by using a parallel development with that of Berdichevsky (1979) 
that the solution obtained here is the first approximation. 

In order to find the stiffness martix, Eq. (18) is used in the first approximation of the 
strain energy functional, Eq. (16) which gives 

U* = - [ e T (A- R T E~ 1 R)edA (19) 

2 Ja 


11 



This is the the first approximation for the strain energy. Derivative of the strain energy 
with respect to the global strain measure e results in the conjugate measures of the section 
stress resultants 


<?=[F Mj r =(^) r (20) 

This relation suffices the existence of a relationship between stress resultants and the 
global measure of strains through an elastic law, which we call the first approximation of 
the matrix of elastic stiffness constants , S*, as in Eq. (10) 


Q = S*e 


( 21 ) 


Following the above operation, S* can then be found as 


S' = A-R t E-'R 


( 22 ) 


The matrix representing the first-order warping contribution to the stiffness matrix, 
R T E~ 1 R , is positive definite. It can also be shown that a finite element approximation for 
the first asymptotic approximation S* will be an upper bound on the actual stiffnesses. 


The First Approximation for Isotropic Plates 

When we reduce our equations to the isotropic case it is possible to obtain an analytical 
solution for warping and the stiffness matrix. Using Eq. (16) for the isotropic case (for the 
first approximation for isotropic case, it is not necessary to discretize the warping) gives 
the warping displacements as follows 


W\ = W2 = 0 


U >3 = —2v 


(711 + 722)2:3 + (kh + /C22) 



(23) 


where v is Poisson’s ratio. The the elastic stiffness constants can be expressed in terms of 
1 /, Young’s modulus E, and the shear modulus G. For the first approximation, the matrix 
of elastic constants is found to be 
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Eft 3 T . 

12(1 —v 2 ) ^ 
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_J?ft 3 
12 (l-v i ) ■ 


( 24 ) 


It ran be seen that for the first approximation our results are exactly the same as those from 
classical isotropic plate theory. Since classical plate theory makes use of the plane stress 
assumption we see that the first approximation of our general functional also coincides 
with the plane stress assumption. If we were to neglect the warping starting from the 
beginning of our analysis then the stiffness values would be overestimated. The important 
point here is that classical isotropic plate theory does include warpingl 

As a first step in developing a numerical procedure, a finite element code has been 
written to evaluate the warping and the elastic constants for the isotropic case. Two-noded 
elements were used with C° continuous shape functions. Results which have been obtained 
coincide almost identically with classical theory, and the agreement becomes much better 
as more elements axe taken. An example output from our finite element program is shown 
in Fig. 3 . The distribution of warping through the thickness due to 711 +722 and /Cn +K22, 
respectively, are plotted in Figs. 4 and 5 . 


Elastodynamics of Beams 

In this section, we describe a simplified case for longitudinal dynamics of a beam. We 
begin with the equations of motion and develop a weak form for mixed space-time finite 
elements. Finally, we present some numerical results. 


Linear Rod Elastodynamic Equations 

The equation of motion for a rod is given as 


F' - P + / = 0 


( 25 ) 


where F is the internal axial force, P is the linear momentum, and / is the force applied 
to the rod. All of these quantities are scalars and axe functions of both space and time. 

The linear rod kinematics is defined by a single displacement variable u. The velocity 
of the rod is obtained by differentiating the displacement with respect to time. Then the 
velocity of any point along the rod generator is 
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V = u (26) 

The strain of the rod can be obtained by differentiating the displacement with respect to 
the spatial variable. So that the strain is 


e = u' (27) 

The constitutive laws for the rod, which connect the strain and velocity to axial force 
and linear momentum, can be simply written as 

F — lie P = mv (28) 

The response of the rod can then be obtained by solving these equations simultane- 
ously. Because of the simplicity of the structure, analytical solutions exist for some simple 
loadings. Thus, some inital and two-point boundary value problems can be used a s bench- 
mark cases in which space-time finite element method based upon above equations can be 
compared with some analytical solutions in the literature. 


Weak Formulation for Space-Time Finite Elements 

It is possible to obtain a weak formulation for a rod by just specializing the weak 
form given by Hodges (1990b). However, since for this simple model we have such a 
simple constitutive law, one may find it useful to satisfy the constitutive relationships 
(algebraic equations) strongly. In this way we do not introduce any more unknowns than 
are necessary. Consequently, the following weak form can be obtained 


II ( Su m ) f+ “ (~ sp + fp ) + 

A(x,t ) 



(—P6u + u8P)\ x + <j> 


( FSu - u<5F)| t 


(29) 


where 8u , 8P , and 8F are test functions. It can be seen that one of the important properties 
of this functional is that none of the unknowns are ever differentiated. All the spatial and 
temporal differentiations are performed over the test functions. Weak forms with this 
property have been termed as “the weakest possible form” by Atilgan (1989). However, by 
allowing differentiation only over test functions yields a form such that the field equations 
now govern the test functions; note that they govern the trial functions in the primitive 
weak form. Since one may assign any function as a test function, the Green functions 
of the field equations could be chosen as test functions. The method using this kind of 
weak form has been called the “boundary element method” in the literature. Therefore, 
even though the weak forms can be same, selection of different shape functions can lead 
different solution strategies. This simple example can show that the differences in finite and 
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boundary element techniques axe superficial; both are coming from the same background. 
More details along these fines and weak forms for theoretical mechanics will be found in a 
paper under preperation. 

Applications 

The initial value problem of a cantilevered rod subjected to a suddenly applied load 
(Heaviside step function) at the free end is considered. This is a classical wave propagation 
problem for which the force and linear momentum are discontinuous. This problem was 
also investigated by Iura et al. (1988) by using a different weak form and by Mansur and 
Brebbia (1984) by using the boundary element method. 

Our space-time element is rectangular. With the weak form, Eq. (29), we have chosen 
the following shape functions. For u , F, and P, constants were chosen in the element 
interior. For the boundary, u is a constant but distinct value from the interior on each 
of the space and time boundaries. On the other hand, F and P are represented by Dirac 
delta functions at the element corners. The test functions 8F and 8P are linear in the 
space and time directions, respectively; and 8u is bilinear in space and time directions. 
Our results are shown in Figs. 6-8 for the displacement, force, and linear momentum. 
Notice the discontinuous quantities are predicted accurately. Similar results for a case with 
initial displacement are shown in Figs. 9 — 11. In both cases, the results match the exact 
solution. 


Future Work 

In the near future we will develop the second approximation for plate analysis. When 
applied to the isotropic plate, this will result in a computational method for generating, 
as a check, the so-called “shear-correction factors.” After validation of the code, we will 
then extend it to treat anisotropic, laminated plate problems. 

For the work on beam elastodynamics, we will continue to expand the capability of 
the analysis to deal with periodic excitation, arbitrary beam deformation, and nonlinear 
problems. 
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Appendix 

Consider a line element through the plate thickness (Fig. 12). The principle of virtual 
work for this filament can be written as follows 



rj-i 

6 s Z a dx$ 



sr T zdx j 


(30) 


where Ss is the virtual displacement of an arbitrary point on the normal line ele- 
ment, <5r is a three-dimensional virtual strain, and Z is the three-dimensional stress 
measure conjugate to the strain, arranged in a 6 x 1 column matrix form Z = 
[ Z\\ Z \2 Z'li Z\ 3 Z23 Z33 J T . The tractions on the lateral surfaces of the line 
elements are written as 


£1 


[Z n 



Z 2 = [Z 22 Z\2 Z 23 J 1 


(31) 
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This principle enforces the weak satisfaction of the three-dimensional equilibrium 
equations, and traction-free upper and lower surfaces. This is clearly analogous to the 
beam St. Venant problem, except that this is for a one-dimensional line element through 
the thickness of a plate. 

Let us decompose the displacement field into warping and rigid body displacement 
components 


s 


= W + U + X 3 



(32) 


where w and u axe 3 x 1 column matrices for warping and the displacement of the reference 
surface, respectively. As outlined in the text, it is possible to find the three-dimensional 
strain as 


r = He + hw,3 +Iiw,i +hw ,2 


(33) 


where matrices 7 -f, J 3 , and I 2 are defined in the text. The stress-strain relationship is 
given as Z = DT. Substitution of Eq. (33) into the principle of virtual work, and discretiz- 
ing warping as w = N(x^)W (£ 1 , 22 ) one can obtain the following system of equations 


A. 

B 

D 

I ■ 

B t 

C 

E 

J 

D r 

E T 

H 

K 

.I T 

J T 

K t 

L. 



(34) 


where 


Q=[F M J 


y °~I h 


N T Z a dx 3 


F=|F n F 12 F 22 F 13 F 23 \ t M = [M n M 12 M 22 J t 


and the matrices are defined as 


(35) 


A = 

/ N T iT DI 3 Ndx 3 B = 

/ N T if DI 2 Ndx, 


Jh 

Jh 

C = 

f N T I?DI 2 Ndx 3 D = 

f N T if DI 3 Ndx 


Jh 

Jh 

H = 

f N T ljDI 3 Ndx 3 1= ; 

f N T lfDHdx 3 


Jh J 

h 

J = 

/ N T I^DHdx 3 K= [ 

Jh Jh 

N T if DHdx 3 


E 


= N T I%'DI 3 Ndx 3 


( 36 ) 
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Eq. (34) can be reduced to one matrix equation governing the warping. However, in order 
to solve that equation, one must find a scalar algebraic equation governed by each stress 
resultant (each element of Q ). Since we have not been able to find such an equation, 
we have changed methods. It may be possible to find a solution for Eq. (34) by using a 
different approach, and this possiblity is still open. 
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Figure 3: An example output from finite element code 
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Figure 4: The distribution of war] 



Figure 5: The distribution of war 
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Figure 6: Displacement distribution in space-time domain due to heaviside step function 
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Figure 7: Force distribution in space-time domain due to heaviside step function 
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Figure 9: Displacement distribution in space-time domain due to initial displacement 
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Figure 10: Force distribution in space-time domain due to initial displacement 
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Figure 11: Linear momentum distribution in space-time domain due to initial displacement 






